function [imagen,X,Z,auxiliar]=BeamFormingActivo3(pulso,eco,pitch,vel_prop,retardo,fs,ResX,ResZ,Xi,Xf,Zi,Zf)

%------
%INPUT
%------
%Pulso: Matriz de N sensores x M1 muestras con el pulso enviado.
%Eco: Matriz de N sensores x M2 muestras con el eco registrado
%pitch: Distancia entre sensores, en metros
%efs: Frecuencia de muestreo (en emisión, Hz)
%rfs: Frecuencia de muestreo (en recepción, Hz)
%retardo: tiempo (s) que transcurre entre el comienzo del pulso y que se comienza a registrar el eco
%vel_prop: velocidad de propagación en el medio
%------
%------
%OPCIONALES
%------
%Zi
%Zf
%Xi
%Xf
%ResX: resolución de imagen deseada en X
%ResZ: resolución de imagen deseada en Z
%------

efs=fs;
rfs=fs;

matlabpool size;
poolFlag=0;
if (ans==0)
    matlabpool open
    poolFlag=1;
end
tic

c=vel_prop;
[M1,N]=size(pulso);
eco=[zeros(ceil(retardo*rfs)+M1,size(eco,2));eco]; %OJO con el M1
[M2,N2]=size(eco);
if N2~=N
    disp('Warning: el pulso y el eco deben usar el mismo número de sensores')
    %Esto no es cierto
end
if efs~=rfs
    disp('Warning: El programa no está pensado para tener distinta frecuencia de muestreo en emisión y recepción')
    %Pero es arreglable
end
if nargin<8
    ResX=N;
    ResZ=1024;
end
if nargin<10
    Xi=0;
    Xf=pitch*(N-1);
end
if nargin<12
    Zi=retardo*c/2;
    Zf=Zi+M2*vel_prop/(2*rfs);
end

z1=retardo*c/2; %Debería ser dividido 2
minX=Xi;
maxX=Xf;
minZ=Zi;
maxZ=Zf;
if (minZ<z1)||(maxZ>M2*vel_prop/(2*rfs))
    disp('Warning: usando profundidades (z) por fuera del rango de la señal')
    if (minZ<z1)
        minZ=z1;
    end
    if (maxZ>(M2*vel_prop/(2*rfs) + z1))
        maxZ=(M2*vel_prop/(2*rfs)+z1);
    end
end
DeltaZ=maxZ-minZ;
DeltaX=maxX-minX;
dZ=DeltaZ/ResZ;
dX=DeltaX/ResX;

if (dX<(c/rfs))
    disp('Warning: muestreo espacial por debajo del límite temporal (X)')
    %En el default, esto significa que los sensores están más separados de
    %lo que deberían!
    %Limito el muestreo espacial al límite temporal que tengo
    %dX=c/rfs;
    %ResX=ceil(DeltaX/dX);
end
if (dZ<(c/rfs))
    disp('Warning: muestreo espacial por debajo del límite temporal (Z)')
    %Limito el muestreo espacial al límite temporal que tengo
    dZ=c/rfs;
    ResZ=ceil(DeltaZ/dZ);
end

%Extiendo pulso y analitizo!
zeros_min=ceil(efs*pitch*N/c);
newM1=M1+2*zeros_min;
%newM1=round(2^(ceil(log(newM1)/log(2))));
pulso_ext=zeros(newM1,N);
pulso_ext((zeros_min+1):(zeros_min+M1),:)=hilbert(pulso);
%Calculo fft de pulso extendido
Fpulso_ext=fft(pulso_ext);
Fpulso_ext=Fpulso_ext(1:(newM1/2+1),:);

%Extiendo eco
newM=M2+2*zeros_min;
%newM=round(2^(ceil(log(newM)/log(2))));
eco_ext=zeros(newM,N2);
eco_ext((zeros_min+1):(zeros_min+M2),:)=hilbert(eco);
%Calculo transformada de Fourier
%Feco_ext=fft(eco_ext,newM);
%Feco_ext=Feco_ext(1:(newM/2+1),:);

imagen=zeros(ResX,ResZ);
auxiliar=[];
for i=1:ResX
    i
    tic
    x=minX+(i-1)*dX;
    col_im=zeros(1,ResZ);
    %aux=zeros(1,ResZ,N);
    for j=1:ResZ
        z=minZ + (j-1)* dZ;
        %Calculo respuesta generada por el pulso asumiendo sistema LIT
        [Fpulso_llegada,s1]=propagacionM21([],0,[x-[0:(N-1)]*pitch],z,c,efs,Fpulso_ext);
        [Feco_esperado,s2]=propagacion12M([],s1,[x-[0:(N-1)]*pitch],z,c,efs,Fpulso_llegada);
        %Calculo ventana de tiempo relevante para el eco registrado
        eco_vent=eco_ext(s2+1:(s2+newM1),:);
        eco_esp=ifft(Feco_esperado,newM1);
        %proyeccion=sum(conj(eco_vent).*eco_esp)./sum(abs(eco_esp).^2);
        proyeccion = productoNormalizado(eco_vent,eco_esp);
        %Calculo imagen en el punto en cuestión
        col_im(j)=sum(proyeccion)/N; %Normalizado para que el máximo sea 1.
    end
    imagen(i,:)=col_im;
    toc
end

if poolFlag==1
    matlabpool close
end
Z=minZ+[0:(ResZ-1)]*dZ;
X=minX+[0:(ResX-1)]*dX;
toc
end


function [Gain,retardo]=respuesta_ij(i,j,x,z,c,d)
%i ó j pueden ser un vector (exclusivamente): fila!
%d1=sqrt((x-d*i)^2+z^2);
%d2=sqrt((x-d*j)^2+z^2);
%retardo=(d1+d2)/c;
%Gain=(1/d1 * 1/d2)';
[G1,r1]=prop((i-1)*d,0,x,z,c);
[G2,r2]=prop((j-1)*d,0,x,z,c);
Gain=G1.*G2;
retardo=r1+r2;
end